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Abstract 

We study the dynamic critical behavior of two algorithms: the Swendsen- 
Wang algorithm for the two-dimensional Potts model with q = 2, 3, 4 and a 
Swendsen-Wang-type algorithm for the two-dimensional symmetric Ashkin- 
Teller model on the self-dual curve. We find that the Li-Sokal bound on the 
autocorrelation time T\at,s ^ const x Ch is almost, but not quite sharp. The 
ratio T\nt,s/CH appears to tend to infinity either as a logarithm or as a small 
power (0.05 ^ p ^ 0.12). We also show that the exponential autocorrelation 
time Texp,£: is proportional to the integrated autocorrelation time Tint,£:- 
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1 Introduction 



Monte Carlo (MC) simulations 0, ^ have become a standard and powerful 



tool for gaining new insights into statistical-mechanical systems and lattice field the- 
ories. However, their practical success is severely limitated by critical slowing-down: 
the autocorrelation time r, which roughly measures the MC time needed to produce 
a "statistically independent" configuration, diverges near a critical point. More pre- 
cisely, for a finite system of linear size L at criticality, we expect a behavior r 
for large L. The power z is a dynamic critical exponent and it depends on both the 
system and the algorithm. 

In single-site MC algorithms (such as single-site Metropolis or heat bath), critical 
slowing-down arises fundamentally from the fact that the updates are local; and have 
a dynamic critical exponent z >2. This is a severe critical slowing-down: the total 
amount of computer work needed to study a (i-dimensional lattice of linear size L 
grows times faster than the naive geometrical factor L'^. To study static and 
dynamic critical behavior we need high-precision data (run lengths ^ lO^r), so in 
practice, it is very difficult to obtain high-precision data for large lattices with local 
algorithms. The geometrical factor L'^ is unavoidable for the usual MC simulations, 
so the elimination (or the reduction) of the critical slowing down is the only way to 
make MC simulations feasible close to a critical point. 

In some cases, a much better dynamic behavior is obtained by allowing nonlocal 
moves, such as cluster flips. In particular, the Swendsen-Wang (SW) algorithm for the 
g-state ferromagnetic Potts model achieves a significant reduction in z compared 



to the local algorithms: numerical experiments suggest that 0.22 < 2; < 1, where the 
exact value depends on the dimensionality of the lattice d and on q (See Table |l|). 
The most favorable case is the two-dimensional (2D) Ising model (g = 2), for which 
the numerical data suggest z = 0.2186 ± 0.0068 but may also be compatible with 
a logarithmic growth [0, |l|, @| . In other cases, the performance of the SW algorithm 
is less impressive (though still quite good): e.g., z = 0.515 ± 0.006 for the 2D 3-state 
Potts model [El, z = 0.876 ± 0.011 for the 2D 4-state Potts model pTl, |2l, and 



z ^ 1 for the 4D Ising model p6| , PQ] . Clearly, we would like to understand why this 
algorithm works so well in some cases and not in others. We hope in this way to 
obtain new insights into the dynamics of non-local MC algorithms, with the ultimate 
aim of devising new and more efficient algorithms. 

There is at present no adequate theory for predicting the dynamic critical behavior 
of an SW-type algorithm. However, there is one rigorous lower bound on z due to 



Li and Sokal |T^. The autocorrelation times of the standard (multi-cluster) SW 
algorithm for the ferromagnetic g-state Potts model are bounded below by a multiple 
of the specific heat: 

T"mt,A/-' ^mt,£, ^cxp > COUSt X Ch ■ (1-1) 

Here TV is the bond density [c.f. (P3|)] in the SW algorithm, S is the energy, and Ch 
is the specific heat. As a result one has 
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Table 1: Numerical estimates of the dynamical critical exponent Zi^t,£ for Potts 
models. We also show the values of the ratio a/u. The exact values of a/i^ have been 
taken from Ref. H, |13|. 



where a and u are the standard static critical exponents. Thus, the SW algorithm 
for the g-state Potts model cannot completely eliminate the critical slowing-down in 
any model in which the specific heat is divergent at criticality. 

The important question is the following: Is the Li-Sokal bound ( p. . 1|) / ( |1.2|) sharp 
or not? An affirmative answer would imply that we could use the bound to predict the 
dynamic critical exponent (s) z given only the static critical exponents of the system. 
There are three possibilities: 

i) The bound ([0| ) is sharp (i.e., the ratio t/Ch is bounded), so that (|1.2|) is an 
equality. 

ii) The bound is sharp modulo a logarithm (i.e., t/Ch ~ log^L for p > 0). 

iii) The bound is not sharp (i.e., t/Ch ~ for p > 0), so that ([L^) is a strict 
inequality. 

Unfortunately, the empirical situation, even for the simplest cases, is far from clear. 
From Table |l| we see that the Li-Sokal bound is apparently not sharp for the Ising 
model in d > 2. However, the differences Ziat,e — <y/i' are much smaller ind = 2. A first 
look at Table |l| reveals that the Li-Sokal bound for the Ising model could only be sharp 
if the autocorrelation time grows like a logarithm |jl2|, ^ |^. Otherwise, the bound 
( p..2| ) would be non-sharp. The Li-Sokal bound for the 3-state Potts model seems to 
be non-sharp: the difference Zint.f — a/z^ = 0.115 ±0.006 is clearly not consistent with 
zero. Finally, the 4-state Potts model is rather peculiar: the naive fit to the data, 
^int,£ = 0.876 ±0.011 p3[, is smaller than the (exactly known) value of a/z/ = 1. The 
explanation of this paradox is that the true leading term in the specific heat has a 
multiplicative logarithmic correction, Ch ~ L\og~^^'^ L [jl9|, ^ ||, |2^. Indeed, a naive 
power-law fit to the specific heat yields a/u = 0.770 ± 0.008 [^, consistent with the 
bound (|1.2| ). The Li-Sokal bound would be sharp modulo a logarithm if r ~ L log^ L 



3 



with p > —3/2. Thus, in d = 2, q = 2 and g = 4 are candidates for a sharp (perhaps 
modulo a logarithm) Li-Sokal bound; but for g = 3 this bound is apparently non- 
sharp, at least if the numerical estimates are to be trusted. It is important, however, 
to proceed cautiously in interpreting these estimates, which are after all only fits 
to data from a limited range of L values (typically L < 1024) and which may be 
significantly biased by corrections to scaling. In particular, it is extremely difficult to 
distinguish numerically between a logarithm and a small power 

The above empirical situation motivated in part the study of the 2D Ashkin-Teller 
(AT) model using a SW-type algorithm ||2l[]. This model interpolates between the 
Ising and the 4-state Potts models. The study of the AT model proved to be useful 
to obtain new insights on the sharpness of the Li-Sokal bound. 

In Section ^ we will define more rigorously the concept of autocorrelation time. 
In Section |^ we will review the SW algorithm for the g-state ferromagnetic Potts 
model and the Li-Sokal bound for such algorithm. In Section ^ we will introduce 
two different SW-type algorithms for the AT model. In Section ^ we will analyze the 
sharpness of the Li-Sokal bound for both cases. Finally, in Section |^ we will consider 
the dynamic critical behavior of the exponential autocorrelation time. 



2 Autocorrelation times 

The goal of a MC simulation is to generate random samples of the spins a = 
{cTj} distributed according to a probability distribution 11. This measure is usually 
Gibbsian 11 = e~'^^"'''/Z, where Ti. is the reduced Hamiltonian of the system and Z is 
the partition function (which is generally unknown). 

To do this, we invent a transition probability matrix P(cr — » a') that tells us 
the probability of going from a spin configuration o" to a new configuration a'. This 
matrix should be ergodic (i.e., from any spin configuration one can reach any other), 
and should leave 11 invariant X]{o-} n(o")P(o" — > cr') = n(cr'). The (usually easier to 
prove) detailed balance Il{a)P{a — > cr') = Il{a')P{a' ^ cr) is a sufficient, but not a 
necessary condition for stationarity. 

We then simulate the Markov process defined by the transition matrix P, starting 
at some arbitrary initial configuration a^^^ . In this way we generate a random sequence 
of configurations Given the two properties above (e.g., ergodicity 

and stationarity of 11), it is guaranteed that this Markov chain converges to the 
equilibrium distribution 11, irrespective of the initial configuration. 

In every MC simulation two different autocorrelation times can be defined. Given 
an arbitrary observable of the spins A{a), we can define its equilibrium normalized 
autocorrelation function: 

^Actually, if we do not have access to many orders of magnitude in L, a power-law with a small 
power p is consistent with a logarithm: ^ I + p log L + O (p^ log'^ L) . 
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Typically, one expects that this function decays exponentially (e.g., pAA{t) ~ e"!*!/"^). 
This motivates the following definition of the exponential autocorrelation time for the 
observable A 

t . , 

'rexp,A = hm sup — j . — - (2.2) 

4_.oo - log \pAA[t)\ 

The autocorrelation time Tgxp.yi measures the relaxation time of the slowest mode 
coupled to the observable A. The exponential autocorrelation time of the system is 
the relaxation time of the slowest mode of the system 

Texp = sup Texp.A (2.3) 



This number can be infinite for some perfectly legal algorithms! p5[. The integrated 



autocorrelation time of the observable A is also defined in terms of paa'- 

^ oo 

Tint,A = 2 ^^^^^) (2.4) 

t=— oo 

There are two fundamental (and distinct) issues in dynamical MC simulations: 

1) Initialization bias. If the initial configuration cr(o) is not "characteristic" of the 
equilibrium probability distribution 11, the first MC configurations are not distributed 
according to 11; they will introduce a systematic error in our MC estimates. To 
eliminate this bias, we have to discard the first Nq iterations such that when t > Nq 
the system has reached the equilibrium distribution 11. But how large should A^o 
be? One answer is given by the exponential autocorrelation time Tgxp ( |2.3| ). This 
autocorrelation time is a worst-case bound on relaxation to equilibrium. For most 
systems, it suffices to discard the first A^^o ^ 20rcxp iterations. In this case, the 
deviation from equilibrium (in the sense) will be at most e~^° (^ 2 x 10^^) times 
the initial deviation from equilibrium. 

2) Autocorrelation in equilibrium. Even if the system is in equilibrium, the MC con- 
figurations are not statistically independent, but highly correlated. The MC estimate 
of the mean value (A) is given by 



^ N 



hj^Aia"')- (2-5) 



where N is the total number of measurements. The variance of such estimator can 
be written as 

^ar(^) = ]^ E - 3) ^ ^^2r,„t,A , (2.6) 



when ^ Tint.A- Thus, Yai{A) is 2rint,yi times larger than in independent sampling. 
This means we need a good determination of the integrated autocorrelation time (|2.4| ) 
in order to obtain a reliable estimate of var(A). In practice, the condition N ^ Tint,A 
is replaced by ^ lO^Tint.A [0- 



5 



Close to a second-order phase transition both autocorrelation times ( p.2| ) / ( p.4|) 
diverge as a power law0 



rexp,A ~ min[L,^]^-p-^ 

rint,A ~ min[L,^]^-''^ 

The dynamic critical exponents Zcxp,A and -Zmt.A are in general different 
Section H) . 



(2.7a) 
(2.7b) 

(See 



3 Swendsen— Wang algorithm for the Potts model 



The SW algorithm |2g] is defined for the ferromagnetic g-state Potts model. 
This model assigns to each lattice site i a spin variable cTj taking values in the set 
{1, 2, . . . , g}. These spins interact through the reduced Hamiltonian 



n 



Potts 



-■Jij^iS„,,a, - 1) , 
{ij) 



(3.1) 



where the sum runs over all the nearest-neighbor pairs (ij), and the all the coupling 
constants are positive Jij > ^{ij)- The Boltzmann weight of a configuration {a} is 
given by 



Potts 



(W) = in 



-1) 



(3.2) 



(ij) 



where Pij = 1 — e~'^'^ , and the partition function is given hj Z = ^{o-} e~^p°"=. The 
idea behind the SW algorithm |]26| , Q is to decompose the Boltzmann weight by 
introducing new dynamical variables Uij = 0, 1 (living on the bonds of the lattice), 
and to simulate the joint model of spin and bond variables by alternately updating 
one set of variables conditional on the other set. The Boltzmann weight of the joint 
model is ^ 

WjointiW}; {n}) = -^Y[[{l-Pij)Sn,,fi+PijSa,,a,Sn,,,l] ■ (3.3) 

(ij) 

The marginal distribution of ( |3.3| ) with respect to the spin variables reproduces the 
Potts-model Boltzmann weight ( p.2|) . The marginal distribution of ( p.3|) with respect 
to the bond variables is the Fortuin-Kasteleyn [T^, |TT[ random-cluster model with 
parameter q: 



n 



{ij): n^j=l 



n (i-p^^) 



{ij): n,j=0 



„C({n}) 



(3.4) 



Close to first-order phase transitions, we expect that the autocorrelation times behave like 
where c is a constant. This exponential growth is due to the tunneling among the 



different pure phases coexisting at the (first-order) phase transition. 
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where C{{n}) is the number of connected components (including one-site components) 
in the graph whose edges are the bonds with riij = 1. 

We can also consider the conditional probabilities of the joint distribution ( |3.3| ). 
The conditional distribution Pbond = E{-\{a}) of the {n} given the {a} is as follows: 
independently for each bond {ij), one sets riij = when o"i 7^ aj, and sets Uij = 
and 1 with probabilities 1 — Pij and pij when cTj = cxj. The conditional distribution 
Pspin = E{-\{n}) of the {a} given the {n} is as follows: independently for each 
connected cluster, one sets all the spins CTj in that cluster equal to the same value, 
chosen with uniform probability from the set {1, 2, . . . , g}. 

The SW algorithm simulates the joint probability distribution ( p.3|) by alternately 
applying the two conditional distributions just described (i.e., Psw = -Pbond -Pspin)- 

Step 1. Generate a new bond configuration {n} from Pbond- 

Step 2. Generate a new spin configuration {a} from Pgpin- 

The performance of this algorithm has been discussed in Section |1| (See also Ta- 
ble IID. Let us now briefly review the proof of the Li-Sokal bound (|1.1| )/ (|1.2| ) for the 
homogeneous Potts model Pij = p |T^. The strategy of the proof is simple: first, we 
choose two "slow" observables: the bond and energy densities 

and compute their autocorrelation functions at time lags and 1. Then, using some 
general properties of reversible Markov chains, we will deduce lower bounds for the 
autocorrelation times rint,A (for A = M,£) and Texp. These will in turn imply lower 
bounds on the dynamic critical exponents Zint,A and z^^p. 

From ( p.3|) we can compute the expectation values of bond variables riij conditional 
to the spin variables. In particular, we can compute the autocorrelation function of 
the bond occupation M at time lag 1 [O, [2^ 



pj^u{l) = 1 - ^^^^,f^^,^ ^ 1 - 7^ ' (3-6) 
pCh + (1 - p)E Ch 

where the energy density E and the specific heat Ch are defined by: 

E =!(£:), CH = ^var(£:). (3.7) 

At criticality, p Pc> Q and ^ -Ec > 0, so the constant A does not vanish. 

The correlation functions of N under Psw = PbondPspin are the same as under the 
positive-semidefinite self-adjoint operator Pg^r = Pspin Pbond Pspin- This implies |]T7[ 
that we have a spectral representation 

PM^it) = I Al*lrfz/(A) , (3.8) 







with a positive measure dv. It follows that 

PMM{t) > p(l)'*l . (3.9) 
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If we introduce the inequalities ( p.6| ) / (3.9) in the definition of the autocorrelation 



times (|2.2|) / (|2.4|) for the bond occupation we obtain the corresponding Li-Sokal 
bounds dTT]) 0. 



Similar bounds for the energy E are obtained by using the identity E{M\{oY) = p£ 
p]7| , p^ . This equation implies the following relation between the autocorrelation 
functions of the bond occupation and energy: 



> PAfM{t) ■ 



This implies in particular, the following set of equalities and inequalities 



(3.10) 



' cxp 



pnm{'^] 



1 ^ TjntX 

2 - PmAK 



(3.11a) 
(3.11b) 



From (3.11) the Li-Sokal bounds ( |1.1|) for the energy follow (and also for the ex- 
ponential autocorrelation time Tgxp). We can also deduce from Eq. ( |3.11a| ) that if 
pj\fj\f{l) 7^ as L — i> oo, then p2 |. 



(3.12) 



4 Swendsen— Wang— like algorithm for the Ashkin— 
Teller model 

The AT model is a generalization of the Ising model to a four-state model. To 
each site i of the lattice we assign two Ising spins Ci = ±1 and Xj = ±1; they interact 
through the Hamiltonian 

n = -jY^ a,aj - J'Y^ r,Tj - (riCTj^Tj (4.1) 

ihj) ihj) 

It can be regarded as two Ising models (with nearest-neighbor couplings J and J') 
interacting via a four-spin coupling K. This model enjoys several symmetries: 1) It 
is symmetric under permutations of {a,T,(yT). This implies that ( [4.1[ ) is symmetric 
under permutations of the couplings {J,J',K). 2) On any bipartite lattice, ( ^TTP is 
symmetric under sublattice fiip of two of the spins (a, r, ar). This means that ( |4.1| ) 
is symmetric under sign fiip of two of the couplings {J, J',K). 

In this paper we are mainly interested in one particular case of the Hamiltonian 
( [4.1|) : the symmetric Ashkin-Teller (SAT) model characterized by J = J'. Although 
we do not know how to solve this model analyticly, we have a fairly good under- 
standing of its phase diagram on a square lattice (See Figure p. The line K = 
(dotted line in Figure |l|) corresponds to two decoupled Ising models; its ferromagnetic 
critical point is denoted by DIs. The line J = K corresponds to the ferromagnetic 
4-state Potts model, and its critical point is denoted by P (J = i^' = ilog3). The 
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line J = —K corresponds to the 4-state Potts antiferromagnet, which in known to be 



non-critical at all temperatures T > ITQ]. The 4-state Potts subspace is depicted 



with dash-dotted lines in Figure |I[ The line J = corresponds to an Ising model 
in the variable err; its ferromagnetic critical point is denoted by Is, and its antiferro- 
magnetic counterpart by APIs. Pinally, when = oo we have another Ising model 
with (T = r, and coupling 2 J. Its ferromagnetic critical point is denoted by Is'. 

There are several critical curves in the SAT model. The most important one is 
the self-dual curve e~^^ = sinh2J. This curve is critical only for K < ;|log3 (solid 
curve in Pigure |lD, and is noncritical for K > ^ log 3 (dashed curve in Pigure |lD. The 
critical part belongs to the universality classes of the conformal field theories with 
conformal charge c = 1. Along this line the critical exponents vary continuously and 
they are known exactly and references therein] . At the point P two more critical 
curves emerge: one goes to the critical point Is and the other one goes to the critical 
point Is' at i^' = cxD. Pinally, there is another critical curve emerging from the critical 
point APIs and pointing towards K = — oo. The exact location of these three curves 
is unknown, although there are good numerical approximations [jT^. These lines are 
believed to belong to the Ising universality class. 



2D Symmetric Ashkin— Teller model 
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Pigure 1: The phase diagram of the square- lattice symmetric Ashkin-Teller model. 
This model is symmetric under J —J, so we only depict the half plane J > 0. 
The points DIs, X2, ZP, and P correspond to the models analyzed in this paper (See 
text). 



We have considered two different cluster algorithms for the AT model. In all cases 
we can assume that 

J, J' > \K\ (4.2) 

This condition is always achievable via the symmetries of the AT Hamiltonian (|4.1| ) 
for a bipartite lattice. 
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Direct algorithm: The idea behind this algorithm is the same as in the standard 
SW algorithm for the Potts model: we decompose the Boltzmann weight by intro- 
ducing new dynamical variables, and we then simulate the joint model by alternately 
updating one set of variables conditional on the other set. As we have two distinct 
sets of Ising spins, we expect to introduce two distinct sets of auxiliary variables: 
rriij^nij = 0, 1 (associated to the spins a, r respectively). The Boltzmann weight of 
the joint model is 



W^joint 



({a,r};{m,n}) = 1 J] [e-2(^+^')5^^^,o5„.^,o + 



^-2K _ g-2J' 



1 _ e-2(J'+i^) _ ^-2{J+K) _^ g-2(J+J')' 



5ai,cTj5n,Tj5m,j,l5n,j,l ■ (4.3) 



Indeed, the Boltzmann weight of the original AT model is the marginal distribution 
of ( [4. 3D with respect to the spin variables. From ( [4.3|) we can read off the conditional 
probabilities needed in the SW procedure (See Ref. for more details): 

Step 1. Update the bond variables {m, n} using the conditional distribution Pbond = 
E(.|{a,r}). 

Step 2. Update the spin variables {a, r} using the conditional distribution Pspin = 
E{-\{m,n}). 

This algorithm reduces on the ferromagnetic 4-state Potts subspace J = > to 
the standard SW algorithm. Wiseman and Domany have introduced essentially 
this same decomposition of the Boltzmann weight using a different derivation as ours. 
They then studied numerically the single-cluster version of this algorithm. Here we 
study the many-cluster ("SW") version. 

Embedding algorithm: The direct algorithm is perfectly legal, but it is somewhat 
complicated to write the computer code for its Step 1 in an efficient way. This 
motivated the introduction of a variant algorithm in which we deal with only one 
kind of spin (a or r) at a time. Consider the Boltzmann weight of a given bond (zj), 
conditional on the {r} configuration (i.e., the r spins are kept fixed): it is 

W^bond(a„ a,- r„ r,) = e-2(^+^-'-^) + [l - e-2(^+^-»-^)] 5.,,.^, . (4.4) 

We can simulate this system of a spins using a standard SW algorithm. The effective 
nearest-neighbor coupling 

J'^ = J + Kr,T, (4.5) 

is no longer translation-invariant, but this does not matter. The key point is that 
the effective coupling is always ferromagnetic, due to the condition ( [4.2|) . An exactly 
analogous argument applies to the {r} spins when the {a} spins are held fixed. The 
embedding algorithm for the AT model has therefore two parts: 
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Step 1: Given the {r} configuration (which we hold fixed), we perform a standard 
SW iteration on the a spins. The probabihty Pij arising in the SW algorithm takes 
the value Pij = 1 — exp[— 2(J + KTiTj)]. 

Step 2: Given the {a} configuration (which we hold fixed), we perform a standard 
SW iteration on the r spins. In this case, pij = 1 — exp[— 2(J' + Kaiaj)]. 



Wiseman and Domany also constructed an embedding version of their single- 
cluster algorithm. Furthermore, they showed that, in the single-cluster context, the 
direct and embedding algorithms define the same dynamics^; only the computer im- 
plementation is different. However, this equivalence does not hold for our many- 
cluster algorithm. In the direct algorithm we have independent clusters of a spins 
and r spins that could be fiipped simultaneously. In the embedding algorithm we 
have at each step only one of the two types of clusters. 

The embedding algorithm, due to its simplicity, is the one used in our MC study of 
the square-lattice SAT model. We expect that it lies in the same dynamic universality 
class as the direct algorithm, on the grounds that one SW hit of {a} followed by one 
SW hit of {r} should be roughly equivalent to one joint hit of {cr, r}. Of course, we 
do not expect the autocorrelation times for the two algorithms to be equal, but we 
do expect them to be asymptotically proportional as the critical point is approached. 

The embedding algorithm does not reduce to the standard SW algorithm on the 
4-state Potts subspace. However, we expect that they do belong to the same dynamic 
universality class. This fact has been numerically verified by comparing the dynamic 
data for the 4-state Potts model with the embedding algorithm to those quoted in 
Ref. [0 (which correspond to the standard SW algorithm). We see that the ratio 
of the two autocorrelation times is more or less constant within statistical errors, 
and that there is no systematic trend as L grows. Thus, we conclude that they are 



proportional in the limit L — > oo [21 



direct 

1.516 ±0.035 (4.6) 



''"int,£: 



—embedded 
'int,^ 



From Table ^ we see that the performance of the embedding algorithm is reasonably 
good on the self-dual curve of the SAT model compared to the local algorithms. 
The direct algorithm for the AT model satisfies the Li-Sokal bound (|1 . 1| ) . The 



proof is a straightforward generalization of the one given in Section ^ for the Potts case 
pT| , Appendix A]. As we expect that both the direct and the embedding algorithms 
belong to the same universality class, we expect that there is a Li-Sokal bound also 
for the embedding algorithm. In Table |^ we show the dynamic critical exponents 
and the values of the ratio a/v for several points on the self-dual curve of the SAT 
model. We have studied two points on the SAT self-dual curve between the critical 
Ising and 4-state Potts models; they are denoted by X2 and ZF (See Figure |l|). X2 
= (J ^ 0.344132, ^ 0.147920) is one of the points considered in p9[, and ZF = 



More precisely, this equivalence holds when the embedding algorithm is defined by making a 
random choice of Step 1 or Step 2 at each iteration. 
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Model 






X2 


0.477 ±0.028 |21| 


^ 0.4183 


ZF 


0.733 ±0.014 |1 


2/3 ^ 0.666667 



Table 2: Numerical estimates of the dynamical critical exponent Zint,£ for two points 
on the self-dual curve of the square-lattice symmetric Ashkin-Teller model. These 
points interpolate between the Ising and 4-state Potts critical models (See Figure |l|). 
We also show the values of the ratio a/u. 



(J ~ 0.302923, K ^ 0.220343) corresponds to a model studied by Zamolodchikov 



and Fateev i30|l. In both cases the Li-Sokal bound is satisfied and the difference 



Zint,£ — a/u is not very large. 

5 Testing the sharpness of the Li— Sokal bound 

The most straightforward (and naive) method to study the sharpness of the Li- 
Sokal bound ( |1 . 1| ) / ( p^ ) is to fit the autocorrelation time to a power-law Tint = 
and compare the dynamic critical exponent 2;int to the ratio a/u. This procedure 
has a clear disadvantage: in some models the leading term of the specific heat is not 
a pure power law (e.g., in the 2D g = 2,4 Potts models). And there is no reason 
why this should not happen for the autocorrelation times as well. This motivates 
the search for a different criterium: The Li-Sokal bound ( |1.1| ) tells us that the ratio 
'rint,£ /Ch > const., so we may try to fit it to different Ansatze.0 In this section we will 
consider only the energy autocorrelation times, as this is one of the slowest modes of 
the system. We have considered three scenarios: 

1. The Li-Sokal bound is sharp: Z[^t^£ = ajv. In this case we expect that the ratio 
T~mt,£ I Ch will converge to a constant as L ^ oo. In particular, we have tried 
three different Ansatze compatible with this scenario 



■7"int,£' 



c 



H 



A 

A + BL'^ (5.1) 
A + B/\ogL 



The second (resp. third) Ansatz corresponds to power-law (resp. additive 
logarithmic) corrections to scaling. 



^We have used a standard weighted least-squares method in all the fits presented here. As a 
precaution against corrections to scaling, we impose a lower cutoff L > Lmin on the data points 
admitted in the fit, and we study systematically the effects of varying Lmin on the parameter 
estimates and on the value. In general, our preferred fit corresponds to the smallest Lmin for 
which the goodness of fit is reasonable (e.g., the confidence level is > 10-20%), and for which 
subsequent increases in Lmin do not cause the to drop vastly more than one unit per degree of 
freedom. 
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2. The Li-Sokal bound is sharp modulo a logarithm: Zint,e = ct/i^x log^ with p > 0. 
In this case, the ratio Ti^t^e/Cn will diverge like log^L. We have considered two 
Ansatze compatible with this scenario 



Tint,£ 



c 



H 



A + B\ogL 
AlosfL 



(5.2) 



3. The Li-Sokal bound is not sharp: -Zint,^: = ot/i' + p with p > 0. In this case the 
ratio Tint,£;/CH will diverge like a power law. This scenario is represented by the 
Ansatz: 

^ = v4L"/"+P (5.3) 



C 



H 



Our numerical results ||2T|, ^ show that there are only two Ansatze that 

are likely to the six models we have considered here (See Tables H-H). These are 



'T\nt,e 



H 



A + BlogL 
ALP 



with p small 



(5.4) 



Thus, the Li-Sokal bound seems to be either sharp modulo a logarithm or non-sharp 
by a small power 0.05 ^ p ^ 0.12. 



Model 


Ansatz 




Parameters 


Lmin 


x' 


DF 


level 


Ising 


A + B/logL 






192 


1.50 


1 


22% 




AlogL + B 






96 


1.47 


4 


83% 




ALP 


p 


= 0.0593(23) 


96 


1.35 


4 


85% 


q = 3 


AlogL + B 






64 


1.93 


3 


59% 




ALP 


p 


= 0.084(2) 


32 


1.72 


4 


79% 


q = 4 


AlogL + B 






16 


1.99 


5 


85% 




ALP 


p 


= 0.119(11) 


16 


1.06 


5 


96% 


X2 


AlogL + B 






16 


1.43 


4 


84% 




ALP 


p 


= 0.051(9) 


16 


1.28 


4 


86% 


ZF 


AlogL + B 






16 


1.03 


4 


91% 




ALP 


p 


= 0.077(12) 


16 


1.23 


4 


87% 



Table 3: We show the fits to the ratio T\-at,elCH for the different models we consider 
in this paper. For each Ansatz, we give the value of Lmm, the value of the the 
number of degrees of freedom [DF), and the confidence level. For the power- law 
Ansatz ALp we also give the estimate of the power p. Note that for the 4-state 
Potts, X2 and ZF models, the error bar of the ratio Ti^t,s/CH was computed using the 
triangle inequality. Thus, the corresponding values are somewhat overestimated. 



In the Ising model we have performed a high-precision MC simulation on lattices 
ranging from L = 4 to L = 512. In all cases we measured 10^rint,£- iterations [^. In 
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Model 


Ansatz 


Parameters 


r 


2 

X 


Dr 


1 1 

level 


Ising 




^int,£: 


= 0.2265(50) 


192 


0.0017 


1 


97% 




AlogL + B 






192 


1.44 


1 


23% 




A log L + B log L 






96 


1.83 


4 


77% 




/ilj^ log 


V 


= 0.0504(23) 


yo 


i.OO 


4 


oO/o 


q = S 






= 0.515(6) 


1 oo 


U.44 


o 


cno/ 
8070 




' [/± log + X? J 






OZ 


i.04 


4 


oz/o 


q = A 


4 T Zi^+ p 




= 0.876(11) 


32 


2.54 


4 


64%) 




lop-"^/^ r/Zl lop- 4- R) 






1 6 


1 53 




77% 

1 1 /(J 




Li+P log"^/^ L 


P 


= 0.153(28) 


128 


1.30 


2 


52% 


X2 




Zmt,£ 


= 0.477(28) 


128 


0.39 


1 


53% 










128 


0.42 


1 


52% 


ZF 




Zmt,S 


= 0.733(14) 


32 


1.48 


3 


68% 




L2/3(AlogL + 5) 






16 


1.53 


4 


82% 



Table 4: We show the fits to the integrated autocorrelation times for the different 
models we consider in this paper. For each Ansatz, we give the value of L^m, the 
value of the x^? the number of degrees of freedom (DF), and the confidence level. 
For the power-law Ansiitze we also give the estimates of the powers {p or Zmt.f)- 



order to increase the statistics we merged our data to that of Baillie and Coddington 
|jl|, 0. From Table ^ we see that the sharp Ansatz ( |5.1| ) is clearly disfavored: with 
a larger L^m it achieves a poorer confidence level. On the other hand, there is no 
significant difference between the other two Ansatze ( |5.2| ) / ( p^ ). This can be seen 
clearly in Figure |^(a). If we look at Table ^, we see that the pure-power-law Ansatz 
is disfavored, as it is the logarithmic Ansatz proposed by several authors [0, |T], 
Both of them need a larger value of Lmin than in the other two Ansatze. Finally, 
the non-sharp-by-a-power Ansatz (|5.3 ) is now slightly favored over the non-sharp- 
by-a-logarithm Ansatz ( |5.2|) , but the difference is again probably not significant. In 
Figure @(b) we see clearly that both fits are clearly favored over the simple power-law 
and the logarithmic growth. 

In the 3-state Potts model we have performed a MC simulation on lattices ranging 
from L = 4 to L = 1024. In all cases we measured 10^rint,£- iterations From 
Table ^ we see that the power-law Ansatz ( |5.3| ) to the ratio Ti^t,£/CH is slightly favored 
over the logarithmic Ansatz ( |5.2| ). However, if we look at Table ^ we observe that 
the Ansatz L'^^^{A\ogL + B) is favored over the non-sharp-by-a-power Ansatz: the 
former achieves a similar value of the confidence level with a smaller value of Lmin- 
In this model the finite-size corrections in the specific heat are very important |^ 
and this is possibly the cause that one scenario is favored when we study rint,^:, and 
a different one is favored when we consider the ratio Ti^t^e/Cu- 

In the 4-state Potts model we performed a MC simulation using the embedding 
algorithm for the SAT model on lattices ranging from L = 16 to L = 1024. In all 



14 



C\2 



C\2 
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CD - 
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10 100 
log L 
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Figure 2: Dynamic critical behavior for the Ising model, (a) Fits to the ratio 
Tint,£ /Ch- The solid line shows the fit to the Ansatz AL^; the dashed line shows 
the fit to the Ansatz A + i?logL. (b) Fits to the autocorrelation time Tij^t,£- The 
dot-dashed curve shows the fits to the Ansatz A log^ L + B log L (which is almost 
identical to the curve coming from the Ansatz AL^ log L) . The dashed curve shows 
the Ansatz A + BlogL; and the solid curve, the Ansatz AL^'"*'^. 



cases (but L = 1024), we measured at least lO^Tmt.f iterations; for L = 1024 we could 
only measure 1500rint,£- iterations . If we look at Table |^ we see that the power- law 
fit is favored over the logarithmic fit: for the same value of the latter gets twice 
as much than the former. However, if we look at the fits of Ti^t,£ (Table we 
conclude that the non-sharp-by-a- logarithm Ansatz is favored over the other two (i.e., 
the non-sharp-by-a-power and the pure power-law scenarios). With a smaller value 
of Lmin it achieves the best confidence level of all three. 

Finally, the models X2 and ZF were studied using the embedding algorithm for 
SAT model. The simulations were performed on lattices ranging from L = 16 to 
L = 512. We performed at least lO'^T^^t.e measurements for the ZF model, and 
3 X 10^rint,£- for the X2 model. In Tables we see almost no difference between the 
fits corresponding to the non-sharp-by-a-power p.3| and the non-sharp-by-a-logarithm 
|5.2| Ansatze. 

In all cases we conclude that the Li-Sokal bound ( |1.1D / ( |1.2D is not sharp by a small 
quantity: either a logarithm or a small power p. Furthermore, there is some kind of 
continuity among the values of the power p: it grows from q ^ 0.05 in the Ising model 
to g ~ 0.12 in the 4-state model, irrespective if we go through the 3-state Potts model 
{p ^ 0.08) or through the SAT self-dual curve (p(X2) ^ 0.05 and j9(ZF) ^ 0.08). It 



15 



is extremely difficult to distinguish between these two scenarios with lattices up to 
L = 1024. We need larger lattices L ^ 1024 in order to solve this question. 

It is intriguing why the Li-Sokal bound is so close to be sharp in two dimensions, 
and clearly non-sharp for d > 2. The explanation is unknown; but a reasonable guess 
is that in d > 2 the energy and the bond occupation are not among the slowest modes 
of the system (as in c? = 2). It would be very interesting to characterize such slow 
modes for d > 2. 



6 Exponential autocorrelation times 

So far we have analyzed the behavior of the integrated autocorrelation time. How- 
ever, it is not guaranteed that this behavior coincides with that of the exponential 
autocorrelation time. In general, we expect that Z[^t^A 7^ -^int.A for most observables 



A [^. This issue can be explained using the scaling relation for the autocorrelation 
function (|2l|) 



PAAit;L)^\t\-P^h, 



Texp,A L 



(6.1) 



where pa is an exponent, Ha is a scaling function, and ^ is the correlation length of 
the system. Summing over the time t, we obtain 

nnt,A ~ rj^p^ (6.2a) 

Zmt,A = i^-PA)ZcKp,A (6.2b) 

Thus, both dynamic critical exponents are different unless pa = 0. If this is the case, 
then the scaling relation ( |6.1|) simplifies to 



PAA{t; L) ^ hA 



(6.3) 



where Ha is another scaling function. This Ansatz can be numerically verified by 
plotting p££ as a function of t/Ti^t,£ and testing if the data coming from different 
values of L collapse well onto a single curve. If so, then p^ = 0. 

We have plotted p££ versus t/rint,£: for the Ising model in Figure ^. Similar plots 
for the 3- and 4-state Potts models can be found in Refs. p2|, |21| . In all cases we find 
that pss is close to a pure exponential and that the Ansatz ( |6.3| ) is satisfied within 
statistical errors (usually we find small corrections to scaling for the smaller values of 
L; e.g., in the Ising case the Ansatz ( |6.3|) is satisfied for L > 64). Thus, the conclusion 
is that the exponential and integrated autocorrelation times for the energy are likely 
to have the same critical exponent 

Zmt,£ = Zexp,£ (6.4) 

for the SW algorithm (Potts model) and for the embedding algorithm (SAT model). 

The value of rexp,£- can in principle be obtained by fitting pgg to an exponential. 
However, this is not an easy task as the MC estimates of p^s for different t are highly 
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correlated ||2T|. One should compute the full covariance matrix for these random 



variables. In Ref. a method was proposed to perform this computation. However, 
high-precision data is needed in order to obtain meaningful results: we need at least 
6 X 10'^rint,£ measurements. For all the models considered here we can nevertheless 
obtain crude estimates of the ratio Tint,^/Tcxp,^:, which are slightly smaller than 1. 



q=2 Autocorrelation Function of E 




int.E 



Figure 3: Plot of pseif) versus t/rint,£- for the Ising model. The different symbols 
denote the different lattice sizes: L = 4 (•), L = 8 (+), L = 16 (x), L = 32 (□), 
L = 64 (^), L = 128 (o), L = 256 (*), and L = 512 (A). We have also depicted the 
line corresponding to a pure exponential peeit) = exp{—t/Tiat,e)- 
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